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We investigate the stability of growing vesicles using the formalism of nonequilibrium thermody- 
namics. The vesicles are growing due to the accretion of lipids to the bilayer which forms the vesicle 
membrane. The thermodynamic description is based on the hydrodynamics of a water/lipid mixture 
together with a model of the vesicle as a discontinuous system in the sense of linear nonequilibrium 
thermodynamics. This formulation allows the forces and fluxes relevant to the dynamic stability of 
the vesicle to be identified. The method is used to analyze the stability of a spherical vesicle against 
arbitrary axisymmetric perturbations. It is found that there are generically two critical radii at 
which changes of stability occur. In the case where the perturbation takes the form of a single 
zonal harmonic, only one of these radii is physical and is given by the ratio ILyjL-,, where Lp is 
the hydraulic conductivity and L-y is the Onsager coefficient related to changes in membrane area 
due to lipid accretion. The stability of such perturbations is related to the value of I corresponding 
to the particular zonal harmonic: those with lower I are more unstable than those with higher I. 
Possible extensions of the current work and the need for experimental input are discussed. 

PACS numbers: 82.20.-w, 05.70.Ln, 87.16.D- 



I. INTRODUCTION 

The determination of the lowest energy configuration 
of a vesicle is one of the most widely studied variational 
problems [IHil. The vesicle is modeled as a closed two- 
dimensional surface in three dimensions with a given en- 
ergy The subscript m denotes "membrane", since in 
reality the surface of the vesicle is a membrane in the form 
of a lipid bilayer There are several different models 
which give different expressions for iSm [H- The earliest, 
and most widely studied, is the Canham-Helfrich-Evans 
approach 0^ , for which the energy is given by 



E,^^-J{2H-CofdA, (1) 

where H is the mean curvature of the surface with area, 
A, and k and Co are constants: the so-called bending 
rigidity and spontaneous curvature respectively. The 
variational calculation might typically consist of finding 
the shape of the surface which minimizes E^^ for a given 
vesicle volume V and surface area A. 

The popularity of this approach to determining vesicle 
shape has much to do with the straightforward way it can 
be posed — not requiring any substantial input regarding 
vesicle structure or composition — and the richness of the 
possible shapes which are found 1]. However these rea- 
sons are also partly responsible why the field has been 
slow to develop: going beyond this description is almost 
certainly going to involve more of the physics of vesicles 
and the resulting analysis may not be so elegant. One 
of the most obvious drawbacks of the variational studies 
is that they are static. They give us a snapshot of the 
shape of the vesicle, but do not tell us how the shape 
evolves with time, or the time taken for any new shape 
to come about. 



Several dynamical studies of vesicles have already ap- 
peared in the literature (8l-[l5|. While these preliminary 
investigations have proved useful in initiating research in 
this area, all have been deficient in some way or other. 
For instance, some of them are not truly dynamic, rely- 
ing partially on the results of the static analysis Q , 
while others assumed that vesicle shapes were restricted 
to spheres or axisymmetric ellipsoids [l^ [l5j . In this pa- 
per we describe a systematic approach to analyzing the 
stability of growing vesicles. 

An obvious question is: what dynamics should be im- 
posed on the system? Following previous treatments, 
the work presented here uses the formalism of linear 
nonequilibrium thermodynamics (LNET). This assumes 
that that the vesicles are macroscopic [1^]. This is rea- 
sonable given the size of vesicles, although it is clear that 
in some circumstances fluctuations will be important [itI - 
[2l] | . However it should be noted that similar assumptions 
were made deriving the form of the energy ([T]); it was 
based on an analogy between the rod-like lipids and ne- 
matic liquid crystals, using the methodology introduced 
by Frank 22] as motivation. This is a macroscopic static 
description. The equivalent macroscopic dynamical de- 
scription will involve nematohydrodynamics [23l - |25j . The 
approach taken here can therefore be seen as a natural 
extension of the static description which leads to ([TJ . 

Unlike the static theory, it is necessary to postulate 
a mechanism which takes the system away from equilib- 
rium, albeit slowly so that LNET holds. One mechanism 
could be temperature change, another could be the ac- 
cretion of lipids onto the surface from the environment. 
We choose to model the latter mechanism, although we 
would expect that much of the formalism constructed will 
be more widely applicable. 

As will become apparent, this study is restricted in two 
ways. Firstly, we concentrate on the stability of defor- 
mations (as opposed to the full dynamics) and secondly. 
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for mathematical and presentational simplicity we focus 
on deformations which take a spherical vesicle to an ar- 
bitrary axisymmetric shape. The outline of the paper is 
therefore as follows: in Section [IT] the thermodynamics of 
a mixture of point-like constituents (water) and rod-like 
constituents (lipids) is reviewed with the aim of identify- 
ing the relevant forces and fluxes — and therefore consti- 
tutive relations — in a discontinuous LNET description of 
vesicle dynamics; in Section Hill the formalism required to 
describe the change of shape of the vesicle, is outlined; in 
Section|lV]thesc two aspects are brought together to pro- 
vide a dynamical description of a vesicle growing due to 
accretion which is used to study the stability of a spher- 
ical vesicle to axisymmetric perturbations. We conclude 
in Section |V] with a review of the methodology of our 
approach and on the prospects for future work. There 
are three technical appendices one for each of Sections 
ininll andlVl 



II. THERMODYNAMICS 

In a previous study 15] which was restricted to de- 
formations between spheres and ellipsoids, the thermo- 
dynamics of vesicle growth was presented in a straight- 
forward but minimal way. Here, the aim is to provide a 
more detailed account. Some of the necessary theory — 
LNET, the rheology of nematics and the study of liquid 
crystals — is already present in the literature, though for 
clarity, certain parts are recapitulated (in light of well- 
known texts) whilst the details of specific calculations are 
provided in Appendix |^ 



A. Linear nonequilibrium thermodynamics 

The application of LNET to membrane systems has 
been studied previously [26l - [29j . however the focus has 
been primarily on transport phenomena with the mem- 
brane treated as a single discontinuity separating two 
regions. In contrast, this paper is concerned with the 
membrane itself, and the deformations which occur dur- 
ing the process of growth due to accretion. It is assumed 
from the outset that the bilayer is closed {i.e. a vesi- 
cle) and that the surrounding solution is sufficiently di- 
lute that lipids only attach to the surface of the exist- 
ing bilayer (and do not form other aggregates). Refer- 
ences [30, [sil have previously considered the aggregation 
of amphiphiles (lipids) for which the chemical potential 
of a given species is taken to be a function of the aggre- 
gation number, that is, the number of molecules of the 
same species in the local neighborhood. A similar mech- 
anism is implicitly considered here by assuming that any 
molecular preference to be part of the bilayer, rather than 
part of the solution, is controlled by chemical potential 
gradients. 

The usual LNET approach [l^ is to consider a sys- 
tem sufficiently close to equilibrium that it can be di- 
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FIG. 1: Vesicle system schematic: the system is formed from 
two distinct phases, dilute water-Iipid solution and the Upid 
bilayer, which are partitioned into three regions, the exterior, 
membrane and the interior, labeled I, II and III, respectively. 
Thermodynamic variables in regions I and III are taken to be 
independent of position; there are no diffusion flows, viscous 
flows or chemical potential gradients. Region II, the mem- 
brane, is considered to have reached equilibrium in the sense 
that the molecules are arranged in the usual bilayer configura- 
tion (shown in the exploded section); "tails" pointing inwards 
and long axis orientated along the surface normal. Changes in 
the fluid resulting from transport in and out of the membrane 
are assumed to be confined to very small areas surrounding 
the membrane boundary, these areas are labeled IV and V 
and are taken to be quasi-stationary, that is, state variables 
may vary with position but on the timescale of changes expe- 
rienced in regions I, II and III, they are independent of time. 
The exterior is taken to behave like a large reservoir while, by 
contrast, it is assumed that there is no net exchange of lipids 
between the membrane and the interior. 

vided into very small sub-systems which are effectively 
homogeneous. Though small, these sub-systems are con- 
sidered mesoscopic, that is, still large enough to define 
thermodynamic variables. It is then possible to choose a 
sufficiently large scale on which the variables that char- 
acterize each small sub-system form a continuous field. 
Each sub-system, and so each point in space on the larger 
scale, is taken to obey the Gibbs relation, which can be 
written in the general form 

Tds = du + pAv — {dg}j' ^ , (2) 

where all variables are now functions of position and 
time. Here, T is the temperature, p is the pressure and 
following the literature we use s, the specific entropy, 
given by S/M, where S is the entropy and M is the 
mass. Similarly, u = U /M is the specific internal energy, 
1/ = V/M is the specific volume, and g — u ~ Ts + pu 
is the specific Gibbs energy. The subscripted brackets, 
{. . .}t, p, are used to indicate that both temperature and 
pressure are held constant. Traditionally, the rate of en- 
tropy production is then written as a sum of thermody- 
namic forces and fluxes which are related by constitutive 
equations. 

For many LNET problems the relevant scale is such 
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that the system effectively comprises a small number of 
uniform (independent of position) regions. In these cases, 
gradients of thermodynamic variables between regions 
are taken to be singular. For such discontinuous [l6| 
systems, thermodynamic forces take the form of differ- 
ences (rather than gradients) and thermodynamic fluxes 
become total flows between regions. In this paper vesicles 
are treated in such a way; the membrane is taken to be a 
separate region of high lipid density — arranged in typical 
bilayer configuration — with discontinuous transitions to 
uniform dilute water-lipid solutions on either side. 

B. Growing vesicles 

Consider the isolated system described in Fig. [1] no 
external forces act, no chemical reactions may take place 
and temperature is taken to be constant throughout. The 
boundaries of each region are characterized by the out- 
ward normal, whereby all internal boundaries allow both 
heat and particle transfer but the external system bound- 
ary is adiabatic. Pressure is assumed to be controlled by 
a piston, shown in gray, which allows heat but not parti- 
cle transfer. Choosing a G {I — V} to label the separate 
regions of the system the total mass contained in a region 
is given by 

M«^^A4"=^ /pfedF, (3) 

k fc 

where pk is the partial mass density, k g {/, w} is used to 
label the components (lipid and water respectively) and 
the integral is over volume V of region a. With this in 
place it is possible to introduce the total mass flux of a 
component out of a region 

d M" f 

j Pk {vk - -dA, V a e {I, II, III}, (4) 

where Vk is the partial velocity of component fc, is 
the velocity of the boundary and AA is the area element 
(aligned along the outward normal) of the surface, A", 
containing region a. For external boundaries the usual 
"no-slip" condition applies: v = = v° . However, in a 
departure from (l6j , internal (permeable) boundaries are 
permitted to move. 

Our approach is to consider a separation of timescales 
between the changes which occur to the membrane and 
the dynamics of the surrounding fluid. Both the inte- 
rior and exterior regions are assumed to equilibrate on a 
timescale much smaller than vesicle growth — that is, they 
are taken to be uniform. Physically, such an assumption 
is plausible in the light of molecular simulations which 
indicate that the pressure difference between two regions 
separated by a bilayer are surprisingly large [32|. There 
are undoubtedly far from equilibrium regimes in which 
the vesicle is changing quickly enough for heterogeneous 



pressure differences to arise, but these situations are not 
considered in this paper. Indeed, for the same reason we 
neglect any flow fields which could arise as a result of 
friction with the moving membrane. As a result of these 
assumptions, local mass fluxes at the boundary to both 
regions (interior and exterior) are taken to be indepen- 
dent of position. By contrast, the membrane, however, 
is not uniform. As indicated in Fig. [1] molecules are as- 
sumed to be orientated in a bilayer fashion, so that for 
any shape other than a sphere, the local configuration of 
lipids (e.g. molecular splay) has an angular dependence. 
We make the assumption that such differences do not me- 
chanically affect the flow of mass, either water or lipids, 
into or out of the membrane. Taking partial velocities 
to be in the direction of the outward normal, local mass 
fluxes, given by p^ (vk — v^) are also assumed to be con- 
stant at the boundary to the membrane region. 

The relative configuration of lipids in the membrane is, 
however, still considered important thermodynamically. 
Indeed, for such a simplified description of vesicles — with 
a uniform interior and exterior, driven by mass fluxes 
which do not vary at different points on the membrane — 
it seems reasonable to expect that any dynamical behav- 
ior (or shape change) will involve averaging some thermo- 
dynamic quantity over the membrane. For example, two 
vesicles of different shapes but equivalent average molec- 
ular splay are anticipated to undergo dynamics driven by 
the same total flows between interior, exterior and mem- 
brane regions. With this in mind, confining the details 
to Appendix |^ we find that the total entropy produced 
in the system is given by 

k a=U 

where /i^ is the chemical potential of component k and 
a bar above a variable denotes the average value taken 
over the boundary to a region, in the sense that 

x-'^^J^xdA, (6) 

for some thermodynamic variable x. Immediately it is 
clear that p,'^ = p,k for a G {I, III}. (The average over the 
boundary to region II is addressed below). It should be 
noted that in order to reach the result ([S]), one important 
assumption has been made: that the entropy produced 
due to the re-alignment of molecules as the vesicle grows 
can be neglected. This is plausible in the context of a 
stability analysis of deformations; any such entropy pro- 
duced as a result of a small perturbation in the shape is 
proportional to thermodynamic forces (such as pressure 
and chemical potential gradients) within the membrane. 
These gradients are considered negligible on the scale of 
pressure/chemical potential differences across the mem- 
brane. 

One might now ask: what role is played by the energy 
which arises from lipid interactions in the membrane? In 
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order to answer this, we suppose that the internal en- 
ergy averaged over the boundary to region II (the mem- 
brane) is well defined thermodynamically, in the sense 
that u^^ = M^^(s'\ P^^, {c"}, -i/j^'). Here, u, s and v, are 
defined ear her, Ck = Mk/M is the concentration (of com- 
ponent k) and = /M is introduced as a specific ex- 
tensive variable characterizing orientation dependent — 
in this case, amphiphilic — interactions. Dropping super- 
scripts for simplicity, this implies 

du = Tds — pdD + jlkdck + fitfid-ip, (7) 

k 

where /ij, is the chemical potential of component k and 



With this in place, it is now possible to examine 
how the chemical potential in ([5]) depends on the other 
thermodynamic variables. Assuming the form flk — 
flkiT,p, {ck}, k) we see that 

dflk = i^kdp + li^kdn + {dpkjr, p, K, (H) 

where temperature has been held constant and partial 
specific quantities Uk and jk are defined in the following 
way 

\9p Jt, {e,^,}, K 



du ' 

^^/T, p, {-c,} 



(8) 



Note that for this system T — T and p — p. Assuming 
that the average specific internal energy is homogeneous 
and of first order in the masses of each component implies 
that the average specific Gibbs energy may be written as 



E 

k 



(9) 



that is, a sum of the usual free energy of a system of 
point-like constituents — arising due to the concentrations 
of components — and another term relating to the nematic 
nature of the molecules. As previously mentioned, we use 
the simplest macroscopic model of energy associated with 
the orientation of lipids in a bilayer: that attributed to 
Canham, Evans and Helfrich. Indeed, ignoring concen- 
tration dependent terms (and re- introducing superscripts 
for clarity), comparison with ([1]) leads to the following 
identifications 



f/jp = K and — 



1 



i2H-CordA, (10) 



as K does not scale with the system size {i.e. it is in- 
tensive). In the above the membrane is taken to be of 
constant thickness I , where the factor preceding the inte- 
gral comes from the fact that the Canham- Helfrich- Evans 
model is an integral over energy per unit area js^ . Note 
also that the factor of M comes from the fact that ip is 
a specific, or per unit mass, quantity. Here, in order to 
be consistent with membrane model literature (see for 
example the early sections of [l| or [s^), the previously 
defined "bar averaging" taken over the boundary to the 
membrane has been replaced by an average over a surface 
bisecting the two layers of lipids. A™, the so-called neu- 
tral surface (33| . Furthermore, this approach is extended 
to all thermodynamic variables: an average taken over 
the exterior of the membrane is the same as an average 
taken over the neutral surface. 



ipk 



dp.k 



(13) 



p, {5,^fc} 



For a discontinuous system such as ours it is useful to first 
define the difference notation (for some thermodynamic 
variable x) 



i^, V a e {I, II, HI}. (14) 



Once again following it is now possible to write an 
equivalent expression to (|lip for small finite differences 
between regions: 



Aa, /3 Hk =l'k (Aq, 13 p) +i^k i^a, P 
+ {Aa, p p,k}T, 



(15) 



Restating ([5]) using the above difference notation and 
then substituting for (fT5)) it can be seen that, for small de- 
formations, the rate of entropy produced will have contri- 
butions from pressure differences, energy differences due 
to membrane deformation and chemical potential differ- 
ences. The resultant expression may then be simplified 
by applying a number of assumptions which are valid for 
the type of vesicle system discussed here. Firstly, the 
mass of lipids accreting to the membrane from the in- 
terior is negligible. That is, in contrast to the exterior, 
the interior is not treated as a reservoir. Secondly, the 
membrane thickness, I, is considered small on the scale 
of the system. This permits us to write the volume of 
the membrane, to first order in small parameter as the 
area of the surface which bisects the membrane multi- 
plied by the thickness (i.e. = /A™ -I- 0(Z^)). Finally, 
we assume that both k, the bending rigidity, and Cq the 
spontaneous curvature, remain constant. For such a case, 
the average density of lipids in the membrane and the av- 
erage ratio of lipids between the inner and outer layers 
must be unchanged. Therefore each new unit of mass 
added to the membrane is assumed to increase the area 
of the surface which bisects the membrane by a constant 
factor. The manipulations are left to Appendix 1X1 where 
it is shown that 
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TcTtot =Ap 

+ 7 



dt 



dV 1 r d£:, 



T, p, K 



(16) 



p^^^ is 



where 7 is the surface tension [46[ and Ap = p 
the pressure difference between the exterior and interior. 
Here, surface area A, volume enclosed V and bending en- 
ergy Em are defined in relation to a single mathematical 
surface taken to bisect the two lipid monolayers. 

It is worthwhile noting here that 7 corresponds to 
the "interfacial free energy" in the sense defined in [2ll ]. 
That is, the change in free energy which results from in- 
creasing (decreasing) the surface by adding (removing) 
molecules at constant density — cf. Eq. (jA23p . This con- 
trasts with increasing the surface area by reducing the 
density of a fixed number of molecules: the resultant 
change in free energy from such an approach is referred 
to as the "elastic free energy". Here, as with the ma- 
jority of studies to date, the effects of elastic free energy 
are neglected and the lipid bilayer is assumed to be ef- 
fectively incompressible due to the separation of energy 
scales between stretching and bending energies Q. As 
pointed out in 21], thermal fluctuations are thought to 
be more important to studies of elastic free energy rather 
than the interfacial energy considered here. 

The result (ITB)) allows us to identify the forces and 
fluxes for this nonequilibrium system, and this is dis- 
cussed further in Section III CI In our previous work (l5j 
these were identified largely through physical arguments 
and from the work of Kedem and Katchalsky p^ . [35| . 
We did use simple thermodynamic relations, but only to 
show that the term involving the membrane energy could 
be absorbed into effective forces — an analysis which is 
recapitulated in Section fll CI The thermodynamic anal- 
ysis is this paper is far more extensive, and as such is 
not easy to compare with that presented in [15i] . This 
is especially true of the identification of the various con- 
tributions to the entropy. In [l^ we considered a fluid 
separated into two homogeneous regions and wrote down 
the sum of entropy changes associated with each region. 
However, here we consider the membrane as a separate 
region. With proper consideration given to conservation 
laws, such a sum is equal in size but opposite in sign to 
the entropy produced by the system. As long as care 
is taken to give the correct interpretation to the vari- 
ous contributions, both treatments agree, but with the 
present paper giving a mathematical justification to the 
form of the dynamical equations used in that was 
not present in that previous discussion. 



C. Effective pressure and surface tension 

This paper examines the stability of deformations away 
from a sphere i.e. does a deformation grow or decay? In 



this context it is possible to simplify (ITBl) by appealing 
once again to the rationale used earlier when neglecting 
the entropy contributions from molecular re-alignments 
within the membrane. It is assumed that the rate of 
change in energy due to small membrane deformations 
is a function of only two time-dependent variables: the 
surface area and volume of the membrane. That is, on 
small timescales, changes in the energy of the membrane 
are dominated by the addition of lipids to the surface 
or by changes in the pressure difference across the mem- 
brane. This point will be re-examined in more detail in 
Section [nil For now, taking E„j — En^^V, A) the rate 
of entropy production can be written in terms of an ef- 
fective pressure difference and effective surface tension. 
That is 



d^n: 

dt 



T, p, K 



dV JAdt 

dEm\ fdAi 



(17) 



OA )y \ dt p 



which implies 



where 



dAl 
dt L 

J 1 ^ K 



7cffi-^i> , (18) 



(Ap)^ff = Ap - 



fOE,, 



\ dV 



(19) 



and 



7cff = 7 



fdEn 



V dA 



(20) 



Equation ([181) identifies the rate of entropy production 
for a growing vesicle as a sum of two pairs of forces and 
fluxes: the flow of volume into the vesicle, coupled to a 
modified pressure difference across the membrane; and 
the rate of area increase of the membrane (due to accre- 
tion of lipids) which is coupled to an modified surface 
tension. Finally, the usual linear constitutive relation 
between the fluxes and forces may be invoked. Making 
contact with llSlI the volume flux is written as 



A\dt J J, p 



= Lp (Ap)^jj + L-,7eff, (21) 



where Lp and 



are Onsager coefflcients deflned "per 
unit area" following the convention of the initial pap ers 
detailing the theory of hydraulic conductivity [s^. l35| . 

In summary, we have shown how the growth of vesi- 
cles can be described thermodynamically, and reduced it 
under given conditions to the study of two-dimensional 
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surfaces whose volume changes according to Eq. (|2T]) . 
The rest of the paper is devoted to analyzing the growth 
such surfaces, and especially to their deviations from a 
spherical shape. 

III. MEMBRANE DEFORMATION 

In this Section, the formalism needed to describe sur- 
face deformations is outlined. The aspects of differential 
geometry required are given in many standard textbooks 
(e.g. (Sal), but fortunately they are also used (in part) 
for the variational treatments which have been so preva- 
lent in previous studies of vesicle behavior. The book by 
Ou-Yang et al. gives a good account of the details and 
the reader is referred to it for a discussion of the mathe- 
matical background. For this reason, Q forms the basis 
of the notation used below. 



z 

z 




FIG. 2: Diagram showing the parametrization of axisymmet- 
ric shapes. 



allow us to make contact with ([T} by writing 



A. Surface geometry 

The two-dimensional surface (embedded in three- 
dimensions) which represents the membrane is defined 
by a vector field r = r(u,v), where u and v parametrize 
the surface. The tangent (vector) space associated with 
each point on the surface is then spanned by vectors 
ri = dr/dqi, where / g {1,2}, qi ^ u and (?2 = v- From 
here, the first fundamental form, or metric, is defined as 

gim = ri-r,n, (22) 

where the inverse metric y*^ is defined such that g^-' gjk = 
SI, with 51. the Kronecker delta symbol. Here, g is the 
determinant of gij , given by 

g ^ Is'Ps^'^gimgp,, (23) 

where is an antisymmetric two-dimensional Levi- 
Civita symbol. The determinant is used to define the 
surface area element 

dA = y/gdudv, (24) 

and the unit normal 



H^--g^=L,,. (27) 

For consistency with ^ and the majority of membrane 
related literature, (P7)) is defined here contrary to the 
usual convention of differential geometry, so that the 
mean curvature of a sphere is positive, -ffsphcrc — ^/R- 

With the basics in place we now focus on an explicit 
example of the above formalism: shapes which are in- 
variant under rotation about the 0— axis. Such axisym- 
metric vesicles have been the some of the most studied 
in variational calculations. The chosen parametrization 
is shown in Figure [2] here it is clear that generic surface 
parameters u and v have been replaced by the familiar 
angles 9 and </>, where < < tt is the inclination, and 
< (j) < 2tt is the azimuthal angle. In this case r has the 
form 

r ^ p{e)p + k{e)k, (28) 

where p — cos (pi+ain (jjj and where i, j and k are the unit 
vectors in the x,y and z directions respectively. Tangent 
vectors are now given by 

rg^ p'{e)p + k'ie)k, and = p(0)0, (29) 



ri X r2 



V9 



(25) 



In order to quantify the curvature of a surface it is 
further necessary to define second derivatives ri„i = 
d^r /dqidqm, where the coefficients of the second funda- 
mental form 



where a dash is used as shorthand for the derivative with 
respect to 9 and (p = — sin + cos . From the above 
it is clear the metric gij is diagonal, and using 
and (|25p it is a simple exercise to show that 



Vg = p\ip'f + ik'f 



1/2 



(30) 



rir. 



(26) 



and 
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p'k - k'p 



ip') + ik'Y 



1/2- 



(31) 



Here, the explicit dependence on 9 of functions p and 
k has been dropped for simphcity. The form of second 
derivatives rgg, rg^, r^g and r^^ can also be calculated 
from which it is seen that Lij is diagonal. Inverting gij 
it then follows from (P7|) that 



H = 



k'p" - k"p' 



k' 



{p'f + {k'f 



3/2 



{p'f + {k'f 



1/2 



(32) 



B. Perturbation theory 



The aim of this paper is to obtain stability conditions 
for deformations which take a vesicle from a sphere to 
an axisymmetric shape. This can be achieved by writ- 
ing the shape dependent terms of constitutive relation 
(I2ip as perturbations from a sphere and then comparing 
terms of equivalent order in the small parameter control- 
ling the perturbation. The details of this are discussed 
in Section IIVI however for now, in anticipation, pertur- 
bative expressions are required for A, F, = / HdA 
and ^2 = / H'^dA, i.e. the geometric terms which arise 
in (1211) . Choosing appropriate forms for functions p{9) 
and fc(^), we write 



de sin 61 



(36) 



where 



f2 1 d / . d 



(37) 



is the ^-dependent part of the Laplacian in spherical po- 
lar coordinates. With this in mind, psp may now be 
computed by expanding e in terms of the zonal harmon- 
ics: 



(38) 



1=2 



Here, the zonal harmonics Yi {6) are the usual spherical 
harmonics with ni = 0. The scale or size of the perturba- 
tion, e, is taken as a common factor of coefficients {a;}. 
For our purposes ao = ai = as Yo{9) and Yi{9) cor- 
respond to spherical growth and translation respectively 



resp 



and 



Using the properties of the zonal harmonics [38 



2tt I dOsinO Yi,{e)YiM = ^i^h. (39) 

'0 



L''Yi{e) = -i{i + i)Yi{e), (40) 



p{e) ^ RsmO {I + e{e)) , 



(33) 



it follows that 



and 



k{e) ^ Rcose{i + e{e)) 



(34) 



where R is clearly the radius of the unperturbed sphere 
and perturbation e{9) — considered small on the scale of 
R — defines the resultant axisymmetric shape. 

Consider the surface area A = 27r /J" ^de, this will 
serve as a template for calculating V, and the de- 
tails of which are confined to Appendix [BJ Substituting 
and dSll) into ^ and ^ gives 



A = 2t:R 



I 



dOt 



l + 2e + e^ + -{e'f+0{e 



(35) 

where the explicit 9 dependence of e has been dropped 
and as before e' = de/d9. Here, terms of order greater 
than have been ignored, a choice which will become 
clear in the next section. Integrating the fourth term on 
the right-hand side gives 



A = 4:TtR^ + e^R^ ^ 



1 + + 



+ 0{e% (41) 



As mentioned earlier, similar steps can be taken to write 
perturbative expressions for ^1, ^2 and F, however, in 
contrast to (|4T|) these results are cumbersome as they 
contain terms of third order in e: for example, the sim- 
plest result, V, is given by 



V ^--kR" + e 



1=2 



i?3 - - - ,0 , M (42) 



h=2l2=2 13=2 

where f{li, I2, h) is related to the square of a Wigner 
3-j symbol. The full details — including results for ^1 and 
^2 — can be found in Appendix IB] 
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IV. STABILITY 

We wish to address the question of vesicle stabihty dy- 
namically, and specificany, to determine when a spherical 
vesicle becomes unstable and undergoes a shape change. 
Stability questions of this kind do not require the full dy- 
namics of the system to be constructed and in our case it 
is sufficient that only two variables are time-dependent: 
R{t), the radius of the unperturbed sphere and e(t), 
which gives the magnitude of the perturbation. The pic- 
ture is the following. We study the growth of a vesi- 
cle which is spherical, but deformed by a small amount 
defined by ([55)1 . The geometry of the perturbation — 
specified by the set {a;} — is fixed (time-independent), 
but the size of the perturbation — specified by e{t) — is 
not. We ask if there is a time (and so an R{t)) at which 
e{t) starts to increase with which will signal an insta- 
bility. 

Before carrying out this analysis, the growth mecha- 
nism needs to be specified, that is, the rate at which 
lipids are added to the membrane needs to be quantified. 
The simplest assumption is that lipids attach themselves 
uniformly over the surface at a constant rate A, such that 



— = A^ 
dt 



A{t) = A{Q). 



At 



(43) 



However, as described in [l5|, when using a growth law 
of the above form it cannot be assumed that R{t) is in- 
dependent of £{t) so a more consistent approach must be 
taken by moving to a new variable r[t)^ defined as the 
radius of a sphere with equivalent surface area. Setting 
A = Anr"^ and comparing to (|^T|) gives 



At this stage we briefly note that writing = 2k^2 — 
2kCq^i-\- nC^A/2, and remembering V = V{r,£), implies 
that Em = Em{r,e{r,V)) = Era{A,V). Therefore the 
term {dE-ax/ Ai\T, p. k which arises in (|16p can indeed be 
written in the form (jl7p and so the constitutive relation 
(pri) is justified in the context of a stability analysis. 

In order to calculate (pij) it is first necessary to write 
the partial derivatives {dEjn/dV)^^ and {dE^/dA)y — 
which arise in the effective pressure and effective surface 
tension respectively — in terms of r and e. As in [T5| we 
have 



dE^ 
dV 



dEm 
dV 



fdE„ 



\ de 



and 



dA 



dE„ 



1 

Sirr 
1 

87rr 1 V 9r 



dr 

fdEn 



fdE^ 



V d£ 



(49) 



where it is worthwhile noting that since the volume V has 
no terms of order e, the membrane energy E^ must be 
taken to 0{e^) to ensure that partial derivative (|48|) has 
the contributions of order e which are necessary to per- 
form a stability analysis. The calculation — to first order 
in e — of partial derivatives {dEjn/dV)j^ and {dEjn/dA)y 
is left to Appendix [C] the results are given by (jC9|) and 
(jClOp respectively. In order to calculate the left-hand 
side of (PT|) , the condition can be used to show that 
dr/dt = Ar/2, which, alongside (|47|) gives the result 



1 



\ " 2 



R = r \ 1 - e — > ar 

\ 1=2 



1 + -KI + 1) 



+ 0(e^) 



(44) 

This can be substituted into the previous expressions for 
^1 and V — (jB14|) and p2t respectively — to give the re- 
sults (|Cip and (|C2p while ^2 remains unchanged. For 
clarity, we re- write geometric terms ^1, ^2 and V in the 
following simplified way 



-(-I 



(50) 

Substituting ([50)1 into (PT|) and using the results of Ap- 
pendix [Cl it is possible to find conditions on the growth 
of the vesicle by equating terms of the same order in ep- 
silon. 



and 



Ci=4^r l + e2^f)-f £3^f)+0(£4) , (45) 



6 = 47r 



V = -Trr^ 
3 



(46) 



(47) 



where the time dependence is contained solely in vari- 
ables r{t) and e{t). 



A. Zeroth order: spherical growth 

At zeroth order in e — equivalent to spherical growth- 
the condition which arises is given by 



2 



2k / ^2 



(2) 



Ap--T 3 



1/(2) 



Cor 



(51) 



It can be seen from the definitions (|46l) and (HTj) that 
the terms proportional to S}^^ /V'^'^^ are dependent on the 
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choice of perturbation: something that should not be the 
case for an e = condition. Therefore in order to make 
([STjl independent of perturbation it is necessary to impose 
the condition = 2Lp/r. This is the same identifica- 
tion which arose in (isj and leads to the same equation 
for spherical growth in the linear regime: Eq. (15) of 
Indeed, if lipid accretion is "turned off" by setting 
A = then the equilibrium condition for spherical vesicles 
is also recovered. In our previous paper, the identifica- 
tion = 2Lp/ r was motivated by knowing the spherical 
equilibrium condition a priori^ here, the equilibrium re- 
sult could have been derived independently by using the 
fact that the zeroth order condition must not rely on 
perturbation choice by definition. 



B. First order 



47rr 



6 =47r 



and 



167r 



[2~l{l + l)] 



3„3 



e a' 



327r 



f{l)[l{l + l){2-l{l + l))] 



(55) 



1 - 



167r 



+ {2-l{l + l))] 



3^3 



£ a 



IGtt 



f{l) [1{1 + 1){2^1{1 + 1))] 



(56) 



At first order in e 



(52) 



which, after some manipulation, can be shown to be of 
the form 



(53) 



Here, — 2Lp/ L~f and 



V = -nr-" 



l+'-L^[2-lil+l)] 
+ ^/(0+O(£^) 



(57) 



respectively, where /(/) = /(/, I, I). Using (|45l)-(|47l) and 
substituting into (|53p we write 



a.| = -(^-rcJ(.-r.J^^^5(0, (58) 



where 



and 



20;(^ -H) - 6/^(; + 1)^ 
Co(3Z2(Z + l)2_6/(/-f l) + 8)' 



(59) 



^(3)^(2) 



(54) 



are critical values of r{t) — the radius of a sphere with 
equivalent area — and correspond to critical values for the 
surface area 47rr^^ and 47rr^^ respectively. Whether the 
surface area of the vesicle is greater than, less than or in 
between the two critical values of the surface area, con- 
trols the sign of the right-hand side of (15^ and therefore 
whether a particular perturbation is stable or unstable. 



Single mode perturbations 



9(1) = 



f{l) + 1)^-61(1 + 1)+ I 
{2-l{l + l))' 



(60) 



It is helpful at this stage to introduce dimensionless quan- 
tities: radii are re-scaled by a factor of Co such that 



Cor, 



Corci and — Corc2 whilst time is 



re-scaled to give r = Xt. Inserting this into the stabil- 
ity condition gives rise to a natural choice for re-scaling 
the Onsager coefficient associated with surface growth, 
= L-^kCq/X. Also, noting that g{l) is always positive 
and (0 is always negative for Z > 2 it can be seen that 
for such perturbations — corresponding to a single mode 
of the zonal harmonics — there is only one critical surface 
area, 47rr^j^ . Taking this into account ((58)) becomes 



In order to understand the implications of ([5^ we 
consider a simplified case: perturbations which corre- 
spond to only one mode of the zonal harmonics, that 
is e{9) — eaiYi{9). Substituting into previous results, 
Eqs. dell), dnSl) and (IbTs)) are reduced to 



a; 



de 
d7 



(f - fci) (r + \rc 



27rL, 



5(0. (61) 



In addition, it should be noted that since f{l, I, I) is 
zero for odd values of /, perturbations which correspond 
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to odd I are constant in time for all radii at first order. In 
order to investigate the stability of odd zonal harmonics 
it is necessary to take the analysis presented to next order 
in £ which is left for future work. Indeed, at a heuristic 
level, it might have been anticipated that the stability 
of harmonic perturbations which are asymmetric about 
6 = are determined at an order greater than symmet- 
ric perturbations, as a better "resolution" is needed to 
differentiate between shapes with lower symmetry. 

In order to interpret the single-mode stability condi- 
tion ([5T|) further, it is first necessary to assume some 
typical values of the constants involved: we estimate a 
value for k of 10~^^J and a value for Co of 10''m^^ 
such that Cpr is of order 1 for a lOOnm spherical vesicle. 
Following |l4] we use a value for Lp, the hydraulic per- 
meability, of 7.5 X 10~^^ms~^Pa~^. What is not known 
is a typical value for L-y, the Onsager coefficient linking 
surface tension to the rate of change in the volume of the 
interior. 

It is possible to deduce an estimate for Lj by recalling 
Eq. (|21[) and comparing the relative contributions that 
both pressure and surface tension terms make to the rate 
of change of volume. In order to do this it is neces- 
sary to estimate an order of magnitude for the effective 
pressure difference (Ap)off and effective surface tension 
7eff. For simplicity we ignore the modifications due to 
the membrane energy and drop the subscript effective. 
An estimate of 7 = 10~"^Nm~^ is provided by [i^ though 
estimating the pressure difference is less clear. Following 
[l3 | we ask: what is the pressure difference which main- 
tains mechanical equilibrium in a growing vesicle? Using 
Eq. (15) of [H], 



de/dr 
0.006 



Ar CoK , 



r 



(62) 



it is possible to use an estimate for A along with those 
previously taken for Co, n, 7 and Lp to compute the 
pressure difference needed to maintain equilibrium in a 
vesicle with r = lOOnm. Taking A = lO^^'s^^ (the mid- 
dle of the range proposed in Q) gives an estimate for the 
pressure difference of 0.1-0.01 bar. Going back to (HU, 
we argue that for such processes the term L^7 will be 
neither negligible nor significantly larger than the term 
LpAp. Indeed, for the purposes of an estimate we require 
that the two terms are the same order. Taking a pressure 
difference of 0.1 bar this assumption implies that Lj is of 
the order lO^^'s^^Pa^^. As a check, we may calculate an 
order of magnitude for the re-scaled quantity using 
the values above. This leads to an order of magnitude 
for of one. For completeness, using these estimates 
the dimensionless analogue for the hydraulic conductiv- 
ity, Lp = LptiCg/X is therefore 7.5. 

Using this estimate for Lj it is possible to plot curves 
which show the stability of particular perturbations at 
different radii. Fig.|3]plots de/dr against f, where e = aie 
ifli > 0), for three specific perturbations: the zonal har- 
monics corresponding to ^ = 2, Z = 4 and 1 = 6. First, 




FIG. 3: (Color online) Graph showing the growth of three 
separate perturbations to a spherical vesicle at different radii. 
Three specific perturbations are shown: the zonal harmon- 
ics corresponding to Z = 2 (solid), I = 4 (dashed) and I = 6 
(dotted). As highlighted in the text, order of magnitude es- 
timates for physical constants have been used such that the 
dimensionless quantities L-y and Lp are taken to be 1 and 



7.5 respectively, 
positive. 



In addition, coefficients ai are taken to be 



we note that all curves cross the axis at = 2CoLp/ L^, 
which in the framework of estimates discussed above is 
15. Below this point, all modes are unstable — that is, 
de/dr is positive — with lowest values of I being the most 
unstable i.e. growing at the fastest rate. Above the criti- 
cal point all modes are stable (decay in time) with lowest 
values of I decaying fastest. However, we may note that 
the inverse is true for perturbations defined opposite to 
those discussed, that is a; — > —ai. These perturba- 
tions are stable beneath (with low I modes decaying 
fastest) and unstable above fci (with low I modes grow- 
ing fastest). 

We note that the curves shown in Fig. [3] are asymmet- 
ric about ■ That is, the scale of perturbation growth 
beneath the critical radius is much larger than the scale 
of decay above the critical radius. Similarly, for ai < 0, 
decay beneath is much larger than growth above ■ 
However, as we will discuss in Section |Vl before a more 
comprehensive analysis may be carried out, the estimates 
of the physical parameters used need to be significantly 
improved. As such, a more detailed discussion of the 
implications of this feature is left for future work. 



2. Ellipsoidal perturbations 

We may ask whether (pT|) can be reduced to our pre- 
vious results for (axisymmetric) ellipsoidal deformations 
. In terms of the dimensionless quantities introduced 
above, the stability condition in that paper is given by 



de 
d7 



f,J(5r~4-6)^ 



15 



287r (ci — ai) ' 



(63) 



where the parametrization used to characterize (axisym- 
metric) ellipsoidal perturbations may be written as 
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Pel = -R (1 + aie) sin 0, 
kc\ — R {1 + cis) cos 6. 



(64) 



Here, the notation of Section fill Bl and Fig. [5] has been 
used: the points on the surface which are an angle 9 
from the positive z-axis are a distance poi from the z-axis 
and feci from the x-y plane. Such a parametrization is 
different to the more general approach taken in the main 
part of this paper, however, it is possible to compare the 
two at first order in e. In order to do so, it is necessary 
to move to variable r, the radius of a sphere of equivalent 
area. Using Eq. (17) of [ll] we see that, to first order in 
e, the surface area of an ellipsoid with parametrization 
(|64l) is given by 



^J -\- that the two parameterizations are equivalent 
to first order in e if 



(71) 



which implies 



a2 



(ci - ai) , 



(72) 



Substituting this expression for a2 into (j6ip and using 
the facts that 5(2) = (5/7) ^5/^ and fc^{l = 2) = -5/6, 
the stability condition for ellipsoidal deformations 
can be recovered. 



l + -(2ai+ci)e + 0(£2) 



(65) 



Setting A — Anr'^ gives 



1 



l-|(2ai + ci)e + 0(e2) 



which can be substituted back into ((64| to give 



(66) 



Pol = r 
feci = r 



l + ^(ai-ci)e + 0(e2) 
l-^(ai-ci)e + 0(e2) 



(67) 



cos 9. 



Here, the radial distance to a point on the surface of the 
deformed shape (ellipsoid) is given by 



\/ (Pcif + {kcif = r 



I + -{ci ~ ai) e {3cos^ 9 - 1) + 0{e' 

(68) 

where taking ¥2(9) = \^J^[icos^9 - l) [H] this may 
be re-written as 



\J {Pcif + {kc 



4 n 



3 V 5 



ici~a,)eY2{9)+0{e^) 

m 

With this in place we may turn to the general axisymmet- 
ric parametrization set out earlier in this paper. Noting 
from dill) that i? = r [l + 0{e'^)] Eqs. §^ and JM]) may 
be written 



p = r [1 + e{9) + 0{e^)] sm9, 
k = r [I + e{9) + Oie^)] cos 61, 



(70) 



where the reader is reminded that the function e(9) is 
of order e. From here it follows from the calculation of 



V. CONCLUSIONS AND DISCUSSION 

The purpose of this paper has been two-fold. Firstly, 
to systematically set up the thermodynamic description 
of vesicle growth, and secondly, to analyze the stability 
of deformations: specifically, to determine when spheri- 
cal surfaces are unstable to small axisymmetric pertur- 
bations. 

While thermodynamic descriptions of vesicle growth 
based on LNET have been discussed previously, the form 
of the fluxes and forces were obtained through physical 
arguments, and not derived from a comprehensive anal- 
ysis of the thermodynamics of a membrane bilayer in an 
aqueous environment. The analysis presented in Section 

In] and Appendix [K\ aims to do just this: it is based on the 
concept of a discontinuous system discussed by de Groot 
and Mazur |l^, although considerably elaborated for a 
vesicle system. It is assumed that as lipids are incorpo- 

-rated into the membrane, the entropy produced — due to 
')che realignments of others — can be neglected. However, 

Jthe energy change that accompanies such accretion can- 
not be ignored. Indeed, the assumption of the Canham- 
Helfrich-Evans form — Eq. ([T]) — replaces the complexity 
of the lipid bilayer by an energy defined in terms of the 
geometrical properties of a single surface. 

In order to investigate the implications of this thermo- 
dynamic approach, we focused on the stability of vesi- 
cles growing due to accretion. In the context of such 
a stability analysis, the energy is assumed to rely on 
two time-dependent variables: surface area A{t) and vol- 
ume enclosed V{t). Here, in analogy to a standard two- 
component system partitioned into two regions, the en- 
tropy produced is characterized by an effective pressure 
difference and effective surface tension. The effective 
pressure is just the normal pressure difference modified 
by the term {dEjn/dV)^^, and the effective surface ten- 
sion is the normal surface tension modified by the term 

idE^/dA)y. 

Whilst undoubtedly such a thermodynamic approach 
can be improved upon, the nature of the approximations 
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made are argued for on physical grounds, and specific 
physical justifications have been given where possible. At 
the very least, the nature and extent of the assumptions 
are clearly visible providing a starting point for attempts 
to relax them. 

With the study reduced to the dynamics of a surface of 
area A{t) enclosing a volume V{t) and characterized by 
an energy £'in(0 — given by ([T]) — the problem is geomet- 
rical in nature and differential geometry may be used. In 
general the approach is applicable to any smooth surface 
with no particular symmetry properties, but the analysis 
that we present in Section IIIII is restricted to an axisym- 
metric surface. The reasons for making this choice are 
two-fold: it is common in variational studies of vesicle 
growth, and, for reasons of simplicity. However there is 
no problem in principle in treating a general shape. The 
main difference is that the perturbation would depend 
on the angle (f) as well as 9, and would be spanned by 
the full set of spherical harmonics Yi,m{&,4')- This would 
lead to sums on I (for instance in Eqs. (|4T|) and (|42|) ) 
being replaced by sums on both I and to. 

When carrying out the stability analysis in Section IIVI 
we further specialized to perturbations about a spheri- 
cal vesicle, but once again other unperturbed geometries 
could be chosen. A peculiarity of the perturbation ex- 
pansion which has already been remarked on in is 
that it is necessary to develop the perturbation expan- 
sion to third order in order to find the growth rate of 
e{t) to leading order. This unfortunately makes the cal- 
culation more complicated than might have naively been 
expected. Nevertheless, the present analysis places no 
constraint on the nature of the perturbation, other than 
it is axisymmetric. This is in contrast to our earlier treat- 
ment [15], which assumed that the distorted surface was 
an ellipsoid. The results of the present treatment are 
shown to reduce to those of [3l for the case where the 
perturbation corresponds to the I — 2 zonal harmonic. 

In general it is found that there are two critical radii 
at which changes in stability occur. However, in the case 
when perturbations take the form of a single zonal har- 
monic, only one of these radii is physical: that given 
by 2Lp/L-y, where Lp is the hydraulic conductivity and 
Lj the Onsager coefficient corresponding to changes in 
surface area due to lipid accretion. In order to quan- 
tify our results we made an estimate for using simple 
physical arguments. For a more complete approach this 
phenomenological coefficient would need to be measured 
experimentally. Using order of magnitude estimates we 
see that — under the conditions of accretion proportional 
to surface area — spherical vesicles at radii other than 
the critical radius are always unstable. Perturbations 
corresponding to zonal harmonics of lowest I grow at the 
fastest rate (are the most unstable) with modes corre- 
sponding to larger I becoming increasingly stable. 

There are many ways in which the analysis presented 
here could be taken forward. The most pressing need is 
for further experimental studies against which the predic- 
tions made in this paper can be compared. Experimen- 



tal work regarding vesicle shape changes has previously 
been carried out though the focus has been on transitions 
induced by either temperature [4lj or osmotic pressure 
changes [42'| . A summary of some relevant experimental 



studies was given in 15[; of particular interest is 42 



where initially symmetric vesicles were found to deform 
into oblate shapes and then prolate ones under osmotic 
pressure changes which serve to reduce the enclosed vol- 
ume. 

In relation to these existing experimental studies, it 
is appropriate to scrutinize our choice of driving mecha- 
nism. Here, we assumed that lipids accrete to the surface 
at a constant rate per unit area. Although this leads to 
a growth law of exponential form (j43p , it is worth high- 
lighting that A is taken to be very small (of the order 
10~^s~^ [^). That is, for the range of observed vesicle 
sizes, surface growth will still be very slow. Our choice 
of mechanism was influenced by a number of factors. 
Firstly, experiments of this type have been carried out 
and are documented in Indeed, the suggestion is that 
growth due to accretion is a plausible phenomenon in the 
pre-biotic scenario in which we are interested. Secondly, 
this approach already entails a great deal of mathemat- 
ics; we therefore wanted to implement the simplest choice 
for driving the system away from equilibrium. However, 
growth laws of a more complicated form could be incor- 
porated if necessary. Finally, we were also conscious to 
demonstrate that both the results presented here and the 
extended theoretical background detailed in Section 
are connected to our previous paper [T5| . Furthermore, 
in reality, lipids in solution are likely to form micelles 
and small vesicles. However, the additional effects of 
such micelle-vesicle or vesicle- vesicle adsorption are not 
considered here. With these points in mind, a more ex- 
perimentally accessible approach might be to extend pre- 
vious work so that temperature change — rather than the 
accretion of lipids — is the effect which gives rise to shape 
changes. However, the analysis needed to incorporate 
such a feature presents a further technical challenge and, 
in light of the already lengthy theoretical background, it 
is left for future work. 

In addition to the above, it is clear from experiment 
that non-axisymmetric vesicles are observed (see for ex- 
ample |43|) and so our analysis should be extended to in- 
vestigate transitions from a sphere to an arbitrary shape, 
as opposed to those invariant about the z-axis. Similarly, 
a more general theory based on deformations from an ar- 
bitrary shape (as opposed to a sphere) could also be de- 
veloped. Lastly, more realistic models of the membrane 
such as the Area-Difference Elasticity model ^4^] could 
be used. 

All these studies would benefit from more experimen- 
tal input. However, the foundations we have laid, and 
the methods we have developed in this paper, do allow 
for these more general analyzes to be carried out. We 
expect that they will lead to a more comprehensive un- 
derstanding of the dynamics of vesicle growth in the near 
future. 
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^ / (Jq + ph {v - v^')) • dA = 0, (A5) 



Appendix A: Entropy production 

In this Appendix, the calculation of entropy produc- 
tion is summarized. For further background the reader 
is referred to chapter 15 of 16] from which this analy- 
sis is adapted. As usual, subscript k £ {w, 1} denotes 
the component (water and lipids respectively) and su- 
perscript a e {I — V} denotes the region (see Figure [ij. 
Starting with the total entropy in each region a 



and 



III 

E 



Pk [Vk - V 



dA = - 



=1 



dt 



= 0, (A6) 



respectively. Imposing (|A5p on (|A4I) and using the facts 
that mass fluxes are assumed to be evenly distributed 
across boundaries and velocities are taken to be in the 
normal direction (see Section |ll]) gives 



psdV, 



(Al) 



the rate-of-change of can be written as a sum of en- 
tropy fluxes over the boundary and any entropy produced 
in the bulk 



~dr 



/ ^(^•^^~'-^n'=^ + ^^^kPk{vk-v''- 



— ph (v — V 



■dA 



(A2) 



adV. 



Here, entropy fluxes have been taken in the traditional 
hydrodynamic form (Equation (20), chapter 3 of (Tg} ) 
plus a term Jncm/T which arises due the nematic nature 
of the lipid molecules [l^ . The other symbols have their 
usual meanings , Jq is heat flow, h = u + is the 
enthalpy, is the velocity of the boundary — zero for all 
impermeable boundaries — and v — VkPk defines both 
total and partial velocities. At this stage, it is assumed 
that entropy fluxes associated with the nematic nature of 
the lipids are tangent to the boundary at any point, and 
so make no contribution to the rate of change of entropy: 



dA = 0. 



(A3) 



In the same manner as set out in chapter 15 of [16| — 
though adapted for the system considered in Fig. [1] — 
Eq. (jA3|) can be used to write an expression for the total 
rate of entropy produced in the system 



Ctot 



^ III . 

a—1 ^ ^ k 



(A4) 



— ph (v — V 



dA, 



where integration is now over internal (permeable) 
boundaries only. Again, adapted from [l6j . conservation 
of internal energy and conservation of mass are given by 



CTtot 



1 



= 1 k 



dM^ 
dt 



(A7) 



where a bar above a variable is used to denote "aver- 
age over a boundary" in the sense of (O. (Note that for 
uniform regions I and III, — f^t)- Using (|A6|) to elim- 
inate mass flows out of region I, the exterior, results in 
([5]) which is re- written here using the difference notation 

(HI 



Ctot 



1 

= y E E ° 



dMl 
dt ' 



(A8) 



As outlined in Section |TT] averages over the boundary to 
the membrane are replaced by averages over the neutral 
surface. Furthermore, it is possible to expand chemical 
potential differences in terms of other thermodynamic 
variables. Using ([TS|) gives 



1 

^tot = ^ E E (^fe (^i. " p) + '^fe (^i^ " 



k a=II 



lA - 1 \dAf," 



(A9) 



So the entropy produced has been written as a sum 
of thermodynamic forces — differences of variables across 
discontinuities — and thermodynamic fluxes. We proceed 
by considering each term separately. The first term is 
summed over regions III and II, the interior and the mem- 
brane respectively. Consider first the interior: a uniform 
region, we may write V^^^ = V^^^ (T, p, {M^^^}, k) and 
therefore 



dV 



III 



dt 



f dV 



III 



T, p, n 



*; / T, p, K, {Af™ } 



dMi" 
dt 



V- TiidAfi" 

k 



(AlO) 
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Consider now the term relating to the membrane (a non- 
uniform region): from Q we have 



where, following from identifications pOj) the second term 
of (jA9|) may finally be written as 



k 



Pk{vk-v'')-AA. (All) 



Recalling that mass fluxes per unit area are assumed 
constant, using the definition of D^^ and the fact that 
X]fe Pk^k — 1 (see Appendix II of [16|]) gives 



T 



T, p, K 



1 

T 



dt 



(A16) 



T, p, K 



In analogy to the Gibbs-Duhem relation derived in Ap- 
pendix II of [T^ , the third term of (jA9p may be simplified 
by writing 



k k -^11 

= / • dA - V / VkPkVk ■ dA (A12) 
Jll f. Jll 

> / VkPkVk ■ dA. 
,. Jll 



dt 



Here, we recognize that the second term on the right- 
hand-side is nothing other than the flow of volume asso- 
ciated with mass moving across a permeable boundary, 
therefore the entire right-hand-side may be written as 
{dV^^/dtjr, p, K- Combining this with ((MT|) the first 
term of (IA9I) becomes 



Q=II J T, p, K 



(A13) 



Turning attention to the second term of (jA9p . we no- 
tice that since interactions between lipids are neglected 
outside the membrane we are free to set k = in all 
other regions, therefore only values of the summand for 
a = II need be considered. In a similar fashion to above, 
remembering that ^ — ipM = J^k ''PkPk, we see that 

k '^^ k -^II 

= / ^i;*" • dA - V / 2/jkPkVk ■ dA 
Ju J, Jll 

- [ - E / ^kPkVk ■ dA. 

Ju 9t ^ Jll 



d^" 



(A14) 

Here, in analogy to above, we recognize the right-hand 
side as {d^'^Vd^lj, ^ ^^ where the second term arises 
due to the fact that ip is not conserved: the curvature 
of the membrane can change spontaneously through the 
exchange of lipids between outer and inner monolayers. 
Finally, assuming that the membrane thickness I is small 
on the scale of the vesicle we may write 



Jll 



^dA + 0{f), 



(A15) 



E '^ki^a, P P'kjr, p, K. = 0, 



(A17) 



from which it can be seen that in the limit of dilute so- 
lutions, Q <C Cu, 



{Ail a M^It, p, « = 0, V a e {I, HI}. (A18) 

Furthermore assuming that the interior of the vesicle only 
contributes a negligible flow of lipids — that is, it is not 
considered a reservoir — gives 



III 



2^ -^{Al, a Pi\t, p, k 



dt 



{Ai_ II fJ,i}T, 



(A19) 

Combining the results (|A13|) . (|A16I) and (|A19p it is now 
possible to write (IA9p as 



Tatot = 2^ (Ai, a p) \ —rr > 

a=II 1^ iT,p,K 



dEn 
dt 



-{Al, 11 Pl}T, p, K- 



T, p, K 



dt 



(A20) 



This expression can be further simplified by once again 
taking the membrane to be of constant thickness, I — very 
small on the scale of the vesicle — so that 



dV^ 
dt 



T, p, K 



dA"' 

"dT 



+ 0(1^), (A21) 



and 



dV 



III 



dt 



T, p, K 



f dV 
\~di 



T, p, K 

I (dA 
2 



dt 



(A22) 



0(1% 



T, p, K 



where A™ and are the area of, and volume enclosed 
by, the surface which bisects the membrane. Finally, as 
K and Co are taken to remain constant, it is assumed 
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that, on average, every unit of mass (of lipids) accreted 
into the bilayer increases the area of the central bisecting 
surface by the same factor, that is 



dMP_ 
dt 



(A23) 



T, p, K 



where a is a constant. This is in-line with the usual 
assumption that bilayers are essentially incompressible 
due to the separation of energy scales between stretching 
and bending energies [2l|. However, it is necessary to ac- 
knowledge that in certain circumstances (e.g. highly com- 
pressed bilayers) thermal fluctuations are important [2l| . 
Applying the results, (|MT|) . (|X22t and (|M3l) to (|X20)) 



gives 



Tatot = (Ai, III p) 

fdA" 

+ 7| 



dV 

~dr 



T, p, K 



dt 



T, p, K 



dt 



T, p, K 



(A24) 



where 



7^(2 (P™ + P') - Ip'' + a{Ai, II w}t, p, .) ■ (A25) 



a = -TT / d^ (pe' + k')^TT d9 (p'e - k') , (B4) 
Jq Jo 



where the second step of (jB4[) comes from integration by 
parts; noticing that p(0) = p(7r) — 0. In order to find 6, 
([33)1 and (pi]) can be substituted into (|B3|) giving 



+ e"-iee'y + l{3e\'-{e'f)' + 0{e% (B5) 



e' = -1 



This expression can be easily integrated. Applying the 
boundary conditions e' (0) = e' (tt) = gives 



+ -ee' + ^ [Seh' - (e') 



(B6) 



Using ([33)) and (p4)) to write down expressions for p' and 
k' and then substituting into (|B4p along with the above 
gives 



S,i=nR de sin ( 
'0 



2 + 2e + {e'f + {e'f ie" - e) 



where the following result has been used 



(B7) 



Appendix B: Perturbative expressions 



In Section llll Bl it is stated that geometrical quantities, 
V, ^1 and ^2 are required in terms of i?, the radius of a 
sphere, and perturbation e{9). This Appendix outlines 
how to arrive at the necessary results. 

Consider first 6 = / HdA = 2tt H^dB, where H 
and y/g are given by ([5^ and ([5D|) respectively. Whilst 
it is possible to directly calculate a perturbative form 
for ^1— by substituting ([33]) and jMl) into ([32]) and ^ 
and integrating — here, a change of variable is introduced 
to simplify the manipulations slightly. Motivated by the 
form of (|30)) introduce variables 7?.(0) and 0(0) such that 



/•TT nlT 

/ de cos e {t'f = -3 / dO sin 9 {e'f t" 
Jo Jo 



(B8) 



A similar procedure may now be applied to ^2 = 
/ H^dA — 2tt H^y/gdO. Using the same variable 
change as above 



it can be immediately seen from earlier definitions (jBl 
that the second term simplifies: 



p'=7^cose, and /^' = 7^sine, (Bl) 
from which it follows that 



7^ 



e' = 



(p'r + (fc'r 



k"p' - p"k' 



and 



(B2) 



(B3) 



dO 



k'B' 



de sine 

=0 

[coseigio = - 



/ "I f'='r 



= 2, 



9=0 



(BIO) 



where the last step follows from the definitions of p', TZ 
and boundary conditions e' (0) = e' (tt) — 0. The re- 
maining terms of (jB9p can be calculated in a straight- 
forward way though the lengthy intermediate steps have 
been omitted here. The result is that 
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^2 = TT / d9 sin t 
Jo 



2- {e'cote + e")- {e'y 

+ i(e'cot0 + e")' + 2(e'f (e + e") 
-e{e'cote + e"f + 0{e^) . 

(Bll) 



In the same fashion as shown in Section IIIIBI for the 
surface area, it is possible to re-write (|B7P and (|B11|) in 
terms of the operator — defined in ([37]) — using inte- 
gration by parts. 



a = 7ri? / d6lsin6' 
Jo 



2 + 2e~ eL^e+ -e^L^e 



(B12) 



and 



^2 = TT / d9 sin ( 
'o 



-e(L2e)'-(e')2i2e + 0(e4) 

"(B13) 

These expressions can then be integrated using ([55]) . 
and (|40|) to give the following results: 



a =47ri?-|-e2^a2Z(/ + 1 



;=2 



e^f a/,a,,a,3<j/3(;3 + l) 2Z2(?2 + 1) 



(B14) 



and 



1=2 



^^=in + e'-Y,af\l{l + l)\l{l + l)~2 



-e^^ H a.i,ai,aiAl3{k + l)\l3ih + l)-2 

(B15) 

where, the function /(/i, ^2, ^3) is given by 



/Gi, ^2, ;3) = 27r / desmeYi,{e)yh{o)YiM 
Jo 



{2h + l){2h + l){2h + l) 



h h I3 




47r 
2 



(B16) 



and 



h h h 




(B17) 



is a Wigner 3-j symbol (see, for example. Appendix C.I 
of (45I I) with m- values set to zero. The symbol is zero 
unless the triangle condition, |^i — Z2I < ^3 < ^1+^2, 
holds. Finally, the volume contained by an axisymmetric 
surface is given by 



2 

= -ttR^ / sin ede [1 + 3e + "ie^ + e^] , (B18) 
3 Jo 



which can also be integrated using the properties of the 
zonal harmonics to give Eq. in the main text. 



Appendix C: Partial derivatives 

In order to write down the right-hand side of (|21l) . 
partial derivatives (|^5)) and are needed in terms of 
r, the radius of a sphere with equivalent surface area, and 
e. This Appendix provides the details of the calculation. 

First, after invoking the growth law (|43|) . the unde- 
formed radius R appearing in the expressions for ^ and 
V — (|B14[) and (|42)) respectively — must be eliminated in 
favor of r. Using the relation pi)) gives 



;=2 



e'^ E ai,ai,ai,\h{h + l)\2h{l2 + l) (CI) 



and 



A 3 



1=2 



(C2) 



^1 ,^2 ,^3 
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where ^2 remains unchanged. We may now follow [IE 
and introduce the reduced volume w, such that 



47rr^ 
V = — - — V. 



(C3) 



Noticing that t; is a function of s only, implies that 



and 



dV J ^ 47rr3/3 



(C4) 



Writing E^^ = 2k^2 - 2kCoCi + kC^A/2, the partial 
derivative [dE^/de)^ can be found directly using (j45|), 
(j46| and A = Aur'^ . Substituting into (fCT]) gives 



V dV 



r3y(2) 



y(2) 



r3T/(2) 



(C9)' 



47rr3/3 



where using the notation of P5)) - (|Tf)) 



where the fact that ^j^^^^^'' = has been used. 

(C5) Similarly, {dE^/dr)^ can also be found from (gSD, ^ 
and A — Attt'^. Substituting this result and (|C9|) into 
gives 



= ^ = 2ey(2) + 3^2^(3) _^ (5(^3)^ (C6) 
de 

Substituting JCH into dm and (IC5|) into dH]) gives 



and 



dE^ 
dV 



47rr3 /a 



9£ 



(C7) 



dA 



kCo nCl 



9k 



■2r2F(2) 



y(2) 



-Car i\ 



:(3) 



1/(3) 



(CIO) 



1 f dEn 



'dEn 



dA J y Sirr \ dr 



3V fdE„ 



2A V dV 



(C8) where once again 

^(2)/y(2) ^ _i/3 has been used. 



[1] U. Seifert, Adv. Phys. 46, 13 (1997). [16 
[2] Z.-C. Ou-Yang, J.-X. Liu, and Y.-Z. Xie, Geometric 

Methods m the Elastic Theory of Membranes in Liquid [17 
Crystal Phases (World Scientific Publishing, Singapore, 

1999). [18 
[3] R. Lipowsky, Nature 349, 475 (1991). 

[4] P. L. Luisi, The Emergence of Life (Chapter 10. Cam- [19 

bridge University Press, Cambridge, 2006). [20 
[5] P. B. Canham, J. Theor. Biol. 26, 61 (1970). 

[6] W. Helfrich, Z. Naturforsch. 28, 693 (1973). [21 

[7] E. A. Evans, J. Biophys 14, 923 (1974). [22 

[8] B. Bozic and S. Svetina, Eur. Biophys. J. 33, 565 (2004). [23 

[9] B. Bozic and S. Svetina, Eur. Phys. J. E 24, 79 (2007). [24 

[10] R. V. Sole and J. Maci'a, J. Theor. Biol. 245, 400 (2007). [25 

[11] R. V. Sole, A. Munteanu, C. Rodriguez-Caso, and J. [26 

Maci'a, Phil. Trans. R. Soc B 362, 1821 (2007). [27 
[12] D. Fanelli and A. J. McKane, Phys. Rev. E 78, 051406 

(2008) . [28 
[13] B. Bozic and S. Svetina, Phys. Rev. E 80, 013401 (2009). [29 
[14] D. FaneUi and A. J. McKane, Phys. Rev. E 80, 013402 [30^ 

(2009) . 

[15] R. G. Morris, D. Fanelli, and A. J. McKane, Phys. Rev. [31] 
E 82, 031125 (2010). 



S. R. de Groot and P. Mazur, N on- equilibrium Thermo- 
dynamics (Dover Publications, New York, 1984). 
M. Wortis, M. Jaric, and U. Seifert, J. Mol. Liq. 71, 195 
(1997). 

W. Helfrich and R.-M. Servuss, Nuovo Cimento D 3, 137 
(1984). 

U. Seifert, Z. Phys. B: Condens. Matter 97, 299 (1995). 
V. Heinrich, F. Sevsek, S. Svetina, and B. Zeks, Phys. 
Rev. E 55, 1809 (1997). 

O. Farago and P. Pincus, Eur. Phys. J. E 11, 399 (2003). 

F. C. Frank, Discuss. Faraday Soc. 25, 19 (1958). 

J. L. Ericksen, J. Rheol. 5, 23 (1961). 

F. M. Leslie, Arch. Ration. Mech. Anal 28, 265 (1968). 

H.-W. Huang, Phys. Rev. Lett. 26, 1525 (1971). 

A. J. Staverman, Trans. Faraday Soc. 48, 176 (1952). 

S. R. De Groot, P. Mazur, and A. Michels, Appl. sci. Res. 
15, 261 (1966). 

J. W. Lorimer, J. Membrane Sci. 25, 181 (1985). 

B. Baranowski, J. Membrane Sci. 57, 119 (1990). 

C. Tanford, The Hydrophobic Effect (John Wiley & Sons, 
New York, 1973). 

J. N. IsraelachviU, D. J. Mitchell, and B. W. Ninham, J. 
Chem. Soc, Faraday Trans. 2 72, 1525 (1976). 



18 



[32] O. H. Samuli OUila, H. J. Risselada, M. Louhivuori, E. 
Lindahl, and S. J. Marrink, Phys. Rev. Lett. 102, 078101 
(2009). 

[33] A. G. Petrov and I. Bivas, Progr. Surf. Sci. 16, 398 
(1984). 

[34] O. Kedem and A. Katchalsky, Biochim. Biophys. Acta 
27, 229 (1958). 

[35] O. Kedem and A. Katchalsky, Trans. Faraday Soc. 59, 
1918 (1963). 

[36] D. C. Kay, Schaum's Outline of Theory and Problems of 
Tensor Calculus (McGraw-Hill, New York, 1998). 

[37] S. A. Safran, J. Chem. Phys. 78, 2073 (1983). 

[38] G. B. Arfken, Mathematical Methods for Physicists (El- 
sevier, Amsterdam, 1985). 

[39] D. Marsh, Chem. Phys. Lipids 144, 146 (2006). 

[40] D. W. R. Gruen and J. Wolfe, Biochim. Biophys. Acta 



2, 572 (1982). 

[41] J. Kas and E. Sackmann, Biophys. J. 60, 825 (1991). 
[42] J. Pencer, G. F. White, and F. R. Hallett, Biophys. J. 

81, 2716 (2001). 
[43] K. Berndl, J. Kas, R. Lipowsky, and E. Sackmann, Eu- 

rophys. Lett. 13, 659 (1990). 
[44] L. Miao, U. Seifert, M. Wortis, and H-G. Dobereiner, 

Phys. Rev. E 49, 5389 (1994). 
[45] A. Messiah, Quantum Mechanics Vol. 2 (North-Holland, 

Amsterdam, Netherlands, 1962). 
[46] Note that in a previous paper ^151], surface tension was 

given by the symbol a. Here the alternative convention 

of using 7 is adopted as a represents the rate of entropy 

production. 



